A New Equation for Describing Heat Conduction by Phonons 
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A new equation, rooted in the theory of Brownian motion, is proposed for describing heat con- 
duction by phonons. Though a finite speed of propagation is a built-in feature of the equation, it 
does not give rise to an inauthentic wave front that results from the application of the hyperbolic 
heat equation (of Cattaneo). Even a simplified, analytically tractable version of the equation yields 
results close to those found by solving, through more elaborate means, the equation of phonon ra- 
diative transfer. An explanation is given as to why both Brownian motion and its inverse (radiative 
transfer) provide equally serviceable paradigms for phonon-mediated heat conduction. 

PACS numbers: 44.10.+1, 05.40.-a, 05.60.Cd, 66.70.+f 



Thermal conduction has been treated, until quite re- 
cently, within the framework of the classical parabolic 
heat equation, associated with the name of Fourier, or 
by using the hyperbolic heat equation, frequently named 
after Cattaneo 0, • The failure of Fourier's equation, 
which becomes apparent when one examines heat con- 
duction on small scales 0, EL HI HI j is easily grasped, for 
it implies an infinite speed of heat propagation; more 
disappointing, however, is the failure, under similar cir- 
cumstances, of Cattaneo's equation, since it does assign 
a finite value to the speed of propagation. Recognition of 
the inadequacies of these equations prompted Majumdar 
to introduce an equation that has come to be known 
as the equation of phonon radiative transfer (EPRT) ; it 
models conduction of heat in thin dielectric films as the 
transport of phonons, in analogy with the linear Boltz- 
mann equation (LBE) used for describingthe transport 
of photons and monoenergetic neutrons |7j. Since even 
simple versions of EPRT are time-consuming Q, one is 
tempted to replace it with an easier alternative that is 
more trustworthy than the equations of Fourier and Cat- 
taneo. Recently, Chen has proposed an approach where 
Cattaneo's equation is grafted onto the collisionlcss form 
of EPRT 0,1^. Chen's recipe, to be labeled here as 
the ballistic-diffusive approximation (BDA), appears to 
combine simplicity with reliability, and his results are 
highly encouraging. An unsatisfactory feature of BDA, 
traceable to Chen's employment of Cattaneo's equation, 
is the occurrence of an artificial wave front in the diffu- 
sive component of the internal energy (and therefore in 
the total internal energy as well) . Another manifestation 
of the distortion introduced by Cattaneo's equation is a 
small but significant disagreement between the tempera- 
ture profiles predicted by EPRT and BDA. 

In this Letter we propose a new equation for describing 
phonon-mediated heat conduction, and present results 
obtained by employing a simpler, analytically tractable 
version of the equation. Since we wish to compare our 
results with those reported by other workers 0, |a> El > we 
will treat a plane-parallel conducting medium. We will 
denote the average speed of sound by v, and the mean 
free path of phonons by 7; the symbol r = l/v will be 
called the relaxation time and 1/r will be denoted by a. 



The complete new equation will be presented after we 
have considered a simplified, exactly soluble variant, and 
examined its performance; meanwhile we will use the la- 
bel NHE for the simpler version, which can be written, 
with d y = d/dy and d y 



dydy, 



d t T(x,t) = (l-e- at ) K d xx T(x,t), 



(1) 



where T(x, t) denotes the temperature at time t at the 
point x, and k = vl/3 is thermal diffusivity. One sees 
immediately that Eq. (1) reduces to PHE, d t T = k8 xx T, 
in the long-time limit (e~ at <C 1). We also point out that 
the substitution 



s = t - 



1 



(2) 



reduces Eq. (1) to d s T = nd xx T. 

To begin with the statement of the problem: A slab 
of thickness L is initially at a uniform temperature To; 
at time t = 0, one face (say that at x = 0) is raised to 
a temperature Ti and is maintained at this temperature 
thereafter, the other face (at x = L) being kept at the 
temperature To; we wish to find the temperature T{x,t) 
and the flux q(x, t), for t > and < x < L. Since the 
problem has already been treated in the past 0, [H E| , 
we will not give more details other than spelling out our 
notation: AT = T\ — To, £ = x/L, t* = at, s* = as, 
9= (T - T )/ AT, (j) = q/(CvAT), Kn = l/L; C stands 
for the specific heat per unit volume. 

It will be convenient to append to 9 and <j) an appro- 
priate suffix (P, H, or N, depending on whether the re- 
sult pertains to the parabolic, hyperbolic, or new heat 
equation). The equations used for calculating 9p(£,t*), 
6ii(£,,t*), (f>p(£,t*), </)#(£,£*) have been stated (with 
some misprints) by previous authors 4, 6). The transfor- 
mation shown in Eq. (2) enables us to express 9n(£,, t*) 
and <f>N(£>t*) as follows: 

= 6p(t,**), (3) 
M£,f) = (l-e-*>p(£*s*)- (4) 

Figure I compares the performance of NHE with the 
classical equations; here the nondimcnsional tempera- 
ture #(£, t*) and the nondimensional heat flux <fr(£, t*) are 
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FIG. 1: Comparison of the spatial profiles of £*) and 
(/>(£, t*) calculated by using the equation of phonon radia- 
tive transfer (EPRT) , the parabolic heat equation of Fourier 
(PHE), the hyperbolic heat equation of Cattaneo (HHE) and 
Eq. (1), the new heat equation (NHE). The data for the curve 
labeled EPRT have been taken from Refs. 
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plotted against £ at t* = 1 in a slab for which Kn = 1 . A 
similar comparison (involving Fourier's equation, Catta- 
neo's equaton, BDA and EPRT) has been made by Chen 
[IJ |(| ; the coordinates of his EPRT plots have been read 
off and included in Fig. 1 to facilitate comparison. Chen 
has rescaled his data so as to harmonize the definitions 
pertaining to the different equations; in the absence of 
such rescaling the EPRT plots show jumps at the bound- 
aries (see below). One sees that the temperature pro- 
file predicted by NHE is in excellent agreement with the 
rescaled EPRT plots; the flux predicted by NHE deviates 
visibly from the EPRT plot (for £ > 0.7), but the other 
two plots (PHE and HHE) are in severe disagreement 
with the EPRT plot for nearly all values of £. 

Curves showing the surface heat flux 0(0, t) are plotted 
in Fig. 2, and compared with reconstructions of Chen's 
plots. The top panel refers to a slab that is so thin that 
transport becomes essentially ballistic — a situation that 
can be better handled by using the extended version of 
Eq. (I) to be presented shortly. In the middle and bottom 
panels of Fig. 2, one notes a persistent difference between 
NHE and EPRT; this is to be expected because the EPRT 



FIG. 2: Comparison of the temporal profile of the heat flux 
at £ = calculated by using EPRT, the PHE, HHE and NHE; 
for an explanation of the labels refer to Fig. 1. The EPRT 
data have been read from Refs. 00 



data originate from plots which have not been rescaled. 

Let us pause briefly to note that NHE provides, in all 
cases, a better description than that furnished by Catta- 
neo's equation. This means that the performance of BDA 
can be improved by discarding the Cattaneo equation in 
favor of NHE. We return now to NHE and investigate 
the prospects of using it as a complete equation in it- 
self; for this purpose, we will not consider samples where 
transport is dominated by ballistic behavior. 

One would expect, in view of some previous investiga- 
tions [3|, |^,|8|], temperature jumps to occur at the bound- 
aries of sufficiently thin samples. To cope with temper- 
ature jumps, we must replace the boundary conditions 
9(0, t*) = 1 and 9(1, t*) = with 9{-£ Q ,t*) = 1 and 
8(1 + £i,t*) = 0, respectively. Postponing the issue of 
how the values of the extrapolated endpoints (xq = L£q 
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FIG. 3: Profiles of #(£, t*) (upper panel) <fi(£, t*) (lower panel) 
as predicted by NHE. The open squares represent thepredic- 
tion of EPRT (for t* = 1 and 10) reported in Ref. g. The 
dotted curve, derived from a plot in Ref. |g is the output (for 
t* = 100) of BDA. 



and x\ = L£i) are to be found, we will estimate their val- 
ues by using the results published by Chen [6j ; it is perti- 
nent to recall that, in the steady state, the diffusive part 
of his solution satisfies Marshak's boundary condition Q : 
%o = x\ — 21/3, or equivalently £o = £i = 2Kn/3. 

Let us consider now a slab for which Kn = 0.1, and fo- 
cus attention on the temperature profiles at the instants 
t* = 1, 10, 100. An examination of Chen's EPRT plot for 
t* = 100, which is almost linear, shows that £o and £i are, 
as one would expect from the above remark concerning 
the boundary conditions, close to 0.07 0. By choosing 
£o = 0.07 = £i, we were led to plots that agreed well, for 
t* = 10 and t* = 100, with the EPRT counterparts in 
Chen's work; but for t* — 1, a sufficiently close fit to the 
EPRT values could not be obtained without changing £o 
to 0.03. To substantiate these statements, we draw the 
reader's attention to the upper panel of Fig. 3, where our 
plots are compared with Chen's BDA curve (t* = 100) 
or with his EPRT plots (t* = 1,10). For the sake of 
avoiding clutter, we have not displayed the BDA values 
for £ > 0.4, and we have also chosen not to reproduce 
Chen's EPRT plot (for t* = 100), which is close to the 
NHE plot for all values of £ and to BDA for £ > 0.2. We 
would like to add here that it did not seem worthwhile to 



improve the fit by fine-tuning the inputs for the extrap- 
olated endpoints. Finally, we show, in the lower panel of 
Fig. 3, our plots for c6(£,i*) for three particular instants 
(t* = 1,10,100), using the same values of £o as those 
mentioned above (0.03, 0.07 and 0.07, respectively). The 
plots for t* = 10 and t* = 100 are in exceedingly good 
agreement with Chen's EPRT plots; the shape of the flux 
plot at t* = 1 accords with its EPRT counterpart, but 
its initial amplitude is larger; we believe that this is due, 
at least in part, to inappropriate choices for the extrap- 
olated endpoints. 

It seems fair to conclude, on the strength of the forego- 
ing evidence, that NHE affords an exceptionally simple, 
stand-alone strategy for studying transient heat conduc- 
tion in a system where precise values of the extrapolated 
endpoints are not needed. On the basis of the results 
presented in Fig. 3, it seems safe to suggest that this 
requirement is satisfied when Kn < 1/10 and t* > 10. 
For other situations, where large jumps at one or both 
boundaries become inevitable, Eq. (1) can still supply 
the diffusive contribution to BDA. 

Having illustrated the performance of an analytically 
tractable form of the new heat equation, we now state its 
general form, 



d t T(x,t) = [(1 



-at\2 



y K d x -u e- at d x T(x,t), (5) 



and hasten to add that the equation is new only in 
the context of heat conduction. If T is replaced with 
F (the probability density in coordinate space), and n 
with D (the diffusion coefficient), Eq. (5) can be iden- 
tified with the equation, derived first by Ornstein and 
van Wijk (O&vW) @, for describing the time evolution 
of the probability density of a Brownian particle that 
starts with an initial velocity uq; their equation reduces, 
if uq has a Maxwellian distribution, to an equation of the 
same form as Eq. (1) [lfj. Our reason for working with 
Eq. (1) was the resulting reduction in labor. Referring to 
their equation, O&vW stated that this is the "diffusion" 
equation which holds for all time intervals; and that it 
reduces, in the limit t — * to d t T = —uoe~ at d x T, de- 
scribing purely ballistic behavior, and to the parabolic 
diffusion equation when t — > oo. Following O&vW, we 
point out that Eq. (5) itself is a special case of the equa- 
tion 

d t T{x, t) = \a{t)d xx T{x, i) - b{t)d x T(x, t), (6) 

whose solution is Gaussian with a mean J Q b(ti)dti and 

a variance J * a(tx)dt\. 

It seems fit to call Eq. (5) the exact ballistic-diffusive 
equation, and to remark that, even in the literature on 
Brownian motion, the analogs of Eqs. (1) and (5) have 
appeared only rarely [lfl, |Tl|. We believe, for reasons 
given below, that Eq. (5) will mimic all the essential fea- 
tures of EPRT. 

A great merit of Majumdar's contribution 3] is that 
it provides a concrete basis for pursuing the analogy be- 
tween phonon-mediated heat conduction and particle dif- 



4 



fusion in terms of the relevant equations. In the origi- 
nal formulation of EPRT, the speed of phonons is taken 
to be a constant, but the mean-free time is viewed as 
an w-dependent quantity, where u> (0 < uj < lod) de- 
notes the frequency. However, in the simplified version 
used in Refs.0, 0, E| and under scrutiny here, the in- 
dependence of the mean free time is ignored. With a 
constant mean free time (denoted here by r) , EPRT can 
be converted, through integration over u>, into an equa- 
tion for g(x,fj,,t) = J Q D [I^(x, fi,t)/v]duj that has the 
same form as LBE with isotropic scattering, which may 
be written as [<9 t + fj,d x ]ip(x ) fj,, t) — a{\V — l)ip(x, fi, t), 

where V = f_-.dfi. This implies that the internal en- 
ergy U{x,t) = T'g{x, [i,t) and the heat flux q(x,t) — 
vVg{x, n,t) are the analogs of F(x,t) = Vijj(x, and 
J(x, t) = vV'4'{x, /x, i), respectively. To accommodate the 
concept of temperature, we recall Chen's definition 
T(x, t) = U(x, t)/C, and stress that a different definition 
would invalidate the rest of the argument. Our next task 
is to show that F^(x, t) may be replaced, to all intents 
and purposes, by F^yi(x,i), where the suffixes RT and 
BM serve to specify the results pertaining to radiative 
transfer and Brownian motion, respectively, and it is to 
this task that we turn now. 

Brownian motion, when modeled as an Ornstein- 
Uhlenbeck process ja, 0] , differs radically from the sit- 
uation portrayed by LBE [l^. However, one is usually 
interested, not in the distribution in phase space, but in 
F(x,t) and J(x,t), and there is good reason to believe 
that the deviation between Fbm and Frt is significantly 
smaller than that between the parent transport equa- 
tions, which are really poles apart. To substantiate this 
claim, we examine two situations: unbounded motion, 
and the Milne problem, where the diffusing particles oc- 
cupy the half-space x > 0, and x = is an absorbing 
boundary. 

Consider a stream of particles which are released, at 
x — Xi with the same initial velocity Uo (parallel to the 
x-axis). It has long been known that Fsu,[x,t;Xi), the 
fundamental solution of the diffusion analog of Eq. (5), 



which was found before the equation itself 12], is a Gaus- 
sian, with a time-dependent mean and a time-dependent 
variance. For particles obeying LBE, Fht(x, t; Xi) can be 
approximated, through an artifice called the cumulant 
approximation (CA), by a Gaussian |14| : here we draw 
attention to the fact that, for isotropic scattering with no 
capture, the expression for F^x{x,t;xi) becomes, when 
CA is imposed, identical with Fsu,{x,t;Xi). 

We move next to the Milne problem, which is ex- 
actly soluble for three linear transport equations [l3j : 
The Klein-Kramers equation, the single-relaxation time 
approximation to the Boltzmann equation, and LBE. 
Though each equation describes a different physical sys- 
tem, the respective plots of F(x, t) against x have the 
same overall appearance; indeed, within the lowest-order 
diffusion approximation, all three imply a rectilinear plot 
with the same extrapolated endpoint; higher-order ap- 
proximations, or the exact solutions, do reveal residual 
differences. Given this background, one expects the tem- 
perature profile predicted by Eq. (5) to be close, but not 
identical, to that found by using EPRT. The results pre- 
sented above bear out this expectation. We would also 
like to mention the possibility that our approach might 
be brought still closer to EPRT by making a choice for 
a(t) and b(t) in Eq. (6) different from that which leads 
to Eq. (5). 

All that remains is to discuss the boundary conditions 
to be imposed on the solutions of Eqs. (1), (5) and (6). 
Since there is some disagreement about this matter in 
the field of heat conduction jfj , we will content ourselves 
by drawing attention to a few relevant contributions con- 
cerning Milne's problem [l5l[T^ . [r^ ]: particularly relevant 
in this context is the finding that the extrapolated end- 
point is a time-dependent quantity. 

It does not seem rash to conclude, on the basis of the 
results and arguments presented above, that the new heat 
equation will provide a very convenient basis for model- 
ing phonon-mediated heat conduction at any spatial or 
temporal scale. 
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